Link to notebook
Link to github repo.
library(tidyverse)
library(phyloseq)
library(dada2)
library(Biostrings)
library(DECIPHER)
library(phangorn)
library(readr)
library(seqinr)
library(decontam)
library(ape)
library(vegan)
#library(philr)
library(RColorBrewer)
library(microbiome)
library(DESeq2)
library(compositions);
library(cowplot)
library(plotly)
library(htmlwidgets)
library(withr)
metadata <- read_csv("sample_data.csv")
[36m──[39m [1m[1mColumn specification[1m[22m [36m─────────────────────────────────────────────────────────────────────────────────[39m
cols(
SampleID = [31mcol_character()[39m,
`Year.Trawl#` = [31mcol_character()[39m,
Datecode = [32mcol_double()[39m,
Date = [31mcol_character()[39m,
Month = [32mcol_double()[39m,
Year = [32mcol_double()[39m,
Bayside = [31mcol_character()[39m,
Station = [31mcol_character()[39m
)
Import count table and taxonomy file. I slightly modified otutable.csv in Excel to otutable_mod.csv to remove the quotes around seq names and put NA placehoder as first col name (which was above row names)
# Import Count table. Skip first row of tsv file, which is just some text
Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls
4: In readChar(file, size, TRUE) : truncating string with embedded nuls
5: In readChar(file, size, TRUE) : truncating string with embedded nuls
6: In readChar(file, size, TRUE) : truncating string with embedded nuls
count_table <- read_table2("results/otutable_mod.csv")
Missing column names filled in: 'X1' [1]
[36m──[39m [1m[1mColumn specification[1m[22m [36m─────────────────────────────────────────────────────────────────────────────────────────────────────────────[39m
cols(
.default = col_double(),
X1 = [31mcol_character()[39m
)
[36mℹ[39m Use [38;5;235m[48;5;253m[38;5;235m[48;5;253m`spec()`[48;5;253m[38;5;235m[49m[39m for the full column specifications.
colnames(count_table)[1] <- "SampleID"
# Import taxonomy of ASVs
taxonomy <- read_csv(file="results/tax_sequences_blast_taxonomy.csv")
Missing column names filled in: 'X1' [1]Duplicated column names deduplicated: 'RefSeq_Tax_ID' => 'RefSeq_Tax_ID_1' [18]
[36m──[39m [1m[1mColumn specification[1m[22m [36m─────────────────────────────────────────────────────────────────────────────────────────────────────────────[39m
cols(
X1 = [32mcol_double()[39m,
ASV_ID = [31mcol_character()[39m,
ref_seq_ID = [31mcol_character()[39m,
PID = [32mcol_double()[39m,
alnmt_len = [32mcol_double()[39m,
mismatch = [32mcol_double()[39m,
eval = [32mcol_double()[39m,
bscore = [32mcol_double()[39m,
RefSeq_Tax_ID = [32mcol_double()[39m,
Ref_Seq_title = [31mcol_character()[39m,
superkingdom = [31mcol_character()[39m,
phylum = [31mcol_character()[39m,
class = [31mcol_character()[39m,
order = [31mcol_character()[39m,
family = [31mcol_character()[39m,
genus = [31mcol_character()[39m,
species = [31mcol_character()[39m,
RefSeq_Tax_ID_1 = [32mcol_double()[39m
)
# remove first col of sequential numbers
taxonomy[,1] <- NULL
# filter out sequences with low PID (recommended by Sara)
taxonomy <- filter(taxonomy, PID > 92)
# remove BLAST metadata and just retain taxonomy (necessary for further processing below)
drop.cols <- c(colnames(taxonomy)[2:9],'RefSeq_Tax_ID_1')
taxonomy <- select(taxonomy, -one_of(drop.cols))
# # And specify that the first column of data are rownames
# taxonomy <- column_to_rownames(taxonomy, var = colnames(taxonomy$ASV_ID))
Filtering removed seqs 110, 332 (Gobiosoma ginsburgi and Belone belone) Note for Sara should we consider setting this at 97% which is more robust and still leaves 334 unique ASVs (rather than 379 with the 92% cutoff in the settings above)
Preview datasets
count_table
taxonomy
metadata
I want to use the phyloseq package for some plotting/ statistics, which first requires making phyloseq objects out of each of input data tables-
count_table_matrix <- as.matrix(count_table[,2:392]) # convert count table to matrix, leaving out character column of sample ID
rownames(count_table_matrix) <- count_table$SampleID # add back in Sample IDs as row names
ASV = otu_table(count_table_matrix, taxa_are_rows = FALSE)
taxonomy_matrix <- as.matrix(taxonomy[,2:8])
rownames(taxonomy_matrix) <- taxonomy$ASV_ID
TAX = tax_table(taxonomy_matrix)
META = sample_data(data.frame(metadata, row.names = metadata$`SampleID`))
First check that the inputs are in compatible formats by checking for ASV names with the phyloseq function, taxa_names
head(taxa_names(TAX))
[1] "Seq_1" "Seq_2" "Seq_3" "Seq_4" "Seq_5" "Seq_6"
head(taxa_names(ASV))
[1] "Seq_1" "Seq_2" "Seq_3" "Seq_4" "Seq_5" "Seq_6"
And check sample names were also detected
# Modify taxa names in ASV, which are formatted with the sample ID, underscor, fastq ID. Don't need this fastq ID anymore and want it to match the sample names from metadata
sample_names(ASV) <- sample_names(ASV) %>%
str_replace_all(pattern = "_S[:digit:]+",replacement = "")
head(sample_names(ASV))
[1] "T1PosCon" "T1S10" "T1S11" "T1S1" "T1S2" "T1S3"
head(sample_names(META))
[1] "T1PosCon" "T1S1" "T1S2" "T1S3" "T1S5" "T1S6"
And make the phyloseq object
ps <- phyloseq(ASV, TAX, META)
rarecurve(otu_table(ps), step=50, cex=0.5)
empty rows removed
# save as .eps
setEPS()
postscript("Figures/rarefaction.eps")
rarecurve(otu_table(ps), step=50, cex=0.5)
empty rows removed
dev.off()
quartz_off_screen
2
Most samples look like they were sampled to completion. Be weary of T3S11, T1S2, and maybe T4S5
Check some features of the phyloseq object
rank_names(ps)
[1] "superkingdom" "phylum" "class" "order" "family" "genus" "species"
unique(tax_table(ps)[, "superkingdom"])
Taxonomy Table: [2 taxa by 1 taxonomic ranks]:
superkingdom
Seq_1 "Eukaryota"
Seq_377 NA
unique(tax_table(ps)[, "phylum"])
Taxonomy Table: [3 taxa by 1 taxonomic ranks]:
phylum
Seq_1 "Chordata"
Seq_368 "Arthropoda"
Seq_377 NA
unique(tax_table(ps)[, "class"])
Taxonomy Table: [5 taxa by 1 taxonomic ranks]:
class
Seq_1 "Actinopteri"
Seq_63 "Mammalia"
Seq_362 "Chondrichthyes"
Seq_368 "Insecta"
Seq_377 NA
There are some ASVs with NA as superkingdom, phylum, or class annotation- delete these.
ps <- subset_taxa(ps, !is.na(superkingdom) & !is.na(phylum) & !is.na(class))
unique(tax_table(ps)[, "superkingdom"])
Taxonomy Table: [1 taxa by 1 taxonomic ranks]:
superkingdom
Seq_1 "Eukaryota"
unique(tax_table(ps)[, "phylum"])
Taxonomy Table: [2 taxa by 1 taxonomic ranks]:
phylum
Seq_1 "Chordata"
Seq_368 "Arthropoda"
unique(tax_table(ps)[, "class"])
Taxonomy Table: [4 taxa by 1 taxonomic ranks]:
class
Seq_1 "Actinopteri"
Seq_63 "Mammalia"
Seq_362 "Chondrichthyes"
Seq_368 "Insecta"
nrow(tax_table(ps)) # number of ASVs left
[1] 378
378 ASVs still remain…
Also check class Mammalia, to see if contamination or real:
tax_table(subset_taxa(ps, class == 'Mammalia'))
Taxonomy Table: [8 taxa by 7 taxonomic ranks]:
superkingdom phylum class order family genus species
Seq_63 "Eukaryota" "Chordata" "Mammalia" "Primates" "Hominidae" "Homo" "Homo sapiens"
Seq_88 "Eukaryota" "Chordata" "Mammalia" "Artiodactyla" "Suidae" "Sus" "Sus scrofa"
Seq_157 "Eukaryota" "Chordata" "Mammalia" "Primates" "Hominidae" "Homo" "Homo sapiens"
Seq_343 "Eukaryota" "Chordata" "Mammalia" "Carnivora" "Felidae" "Felis" "Felis catus"
Seq_369 "Eukaryota" "Chordata" "Mammalia" "Artiodactyla" "Bovidae" "Bos" "Bos taurus"
Seq_378 "Eukaryota" "Chordata" "Mammalia" "Primates" "Hominidae" "Homo" "Homo sapiens"
Seq_383 "Eukaryota" "Chordata" "Mammalia" "Primates" "Hominidae" "Homo" "Homo sapiens"
Seq_389 "Eukaryota" "Chordata" "Mammalia" "Primates" "Hominidae" "Homo" "Homo sapiens"
These are human, wild boar, cat (…cat lady), and cattle. All are contamination so delete all Mammalia
ps <- subset_taxa(ps, !class == 'Mammalia')
unique(tax_table(ps)[, "class"])
Taxonomy Table: [3 taxa by 1 taxonomic ranks]:
class
Seq_1 "Actinopteri"
Seq_362 "Chondrichthyes"
Seq_368 "Insecta"
Next check the “Insecta” entries
tax_table(subset_taxa(ps, class == 'Insecta'))
Taxonomy Table: [2 taxa by 7 taxonomic ranks]:
superkingdom phylum class order family genus species
Seq_368 "Eukaryota" "Arthropoda" "Insecta" "Hymenoptera" "Formicidae" "Linepithema" "Linepithema humile"
Seq_380 "Eukaryota" "Arthropoda" "Insecta" "Hymenoptera" "Formicidae" "Linepithema" "Linepithema humile"
The onlly Insecta is Linepithema humile, which are ants so delete these too..
ps <- subset_taxa(ps, !class == 'Insecta')
unique(tax_table(ps)[, "class"])
Taxonomy Table: [2 taxa by 1 taxonomic ranks]:
class
Seq_1 "Actinopteri"
Seq_362 "Chondrichthyes"
Check overall how the phyla are distributed among samples
# First aglomerate the ASVs at the phylum level using the phyloseq function, tax_glom
superkingdomGlommed = tax_glom(ps, "superkingdom")
# and plot
plot_bar(superkingdomGlommed, x = "Sample")
ggsave(filename = "Figures/seqdepth.eps", plot = plot_bar(superkingdomGlommed, x = "Sample"), units = c("in"), width = 9, height = 6, dpi = 300, )# and save
Total sequences reveals certain samples had very low sequencing effort: T1S7, T1S8, T3S11, and, not as bad, T1S2 and T4S5
The rarefaction analysis also showed T1S2 and T4S5 samples were likely not sequenced to completion. Therefore remove these 5 samples from analysis
ps <- subset_samples(ps, !SampleID == "T1S7" & !SampleID == "T1S8" & !SampleID == "T3S11" & !SampleID == "T1S2" & !SampleID == "T4S5")
ps
phyloseq-class experiment-level object
otu_table() OTU Table: [ 368 taxa and 50 samples ]
sample_data() Sample Data: [ 50 samples by 8 sample variables ]
tax_table() Taxonomy Table: [ 368 taxa by 7 taxonomic ranks ]
50 samples remaining with 368 ASVs
Remove Pos Controls (all hits in positive controls are the same family- I assume this is expected)
ps <- subset_samples(ps, !SampleID == "T1PosCon" & !SampleID == "T2PosCon" & !SampleID == "T3PosCon")
ps
phyloseq-class experiment-level object
otu_table() OTU Table: [ 368 taxa and 47 samples ]
sample_data() Sample Data: [ 47 samples by 8 sample variables ]
tax_table() Taxonomy Table: [ 368 taxa by 7 taxonomic ranks ]
And lastly, correct some taxonomy: According to Sara, Engraulis encrasicolus (European anchovy) should be Anchoa mitchilli (Bay anchovy):
tax_table(ps) <- gsub(tax_table(ps), pattern = "Engraulis encrasicolus", replacement = "Anchoa mitchilli")
47 samples remainwith 368 unique ASVs
For plotting, use relative abundances (# of ASV sequences/sum total sequences in sample), calculated easily using microbiome::transform
ps_ra <- microbiome::transform(ps, transform = "compositional")
Export the relative abundance matrix so Sara can have it:
# Extract abundance matrix from the phyloseq object
RelAbun_matrix = as(otu_table(ps_ra), "matrix")
# Coerce to data.frame
RelAbun_dataframe = as.data.frame(RelAbun_matrix)
# Export
write.csv(RelAbun_dataframe,"results/otutable_relabun.csv", row.names = TRUE)
Then aglomerate the ASVs at the family level using the phyloseq function, tax_glom
familyGlommed_RA = tax_glom(ps_ra, "family")
family_barplot <- plot_bar(familyGlommed_RA, x = "Sample", fill = "family")
family_barplot
NOTES for Sara
Agglomerate by species to see if I get the same 38 unique species Sara sees:
speciesGlommed_RA = tax_glom(ps_ra, "species")
speciesGlommed_RA
phyloseq-class experiment-level object
otu_table() OTU Table: [ 42 taxa and 47 samples ]
sample_data() Sample Data: [ 47 samples by 8 sample variables ]
tax_table() Taxonomy Table: [ 42 taxa by 7 taxonomic ranks ]
tax_table(speciesGlommed_RA)
Taxonomy Table: [42 taxa by 7 taxonomic ranks]:
superkingdom phylum class order family genus
Seq_1 "Eukaryota" "Chordata" "Actinopteri" "Atheriniformes" "Atherinopsidae" "Menidia"
Seq_2 "Eukaryota" "Chordata" "Actinopteri" "Clupeiformes" "Clupeidae" "Brevoortia"
Seq_3 "Eukaryota" "Chordata" "Actinopteri" "Clupeiformes" "Engraulidae" "Engraulis"
Seq_4 "Eukaryota" "Chordata" "Actinopteri" "Scombriformes" "Pomatomidae" "Pomatomus"
Seq_5 "Eukaryota" "Chordata" "Actinopteri" "Lutjaniformes" "Lutjanidae" "Lutjanus"
Seq_6 "Eukaryota" "Chordata" "Actinopteri" "Pleuronectiformes" "Paralichthyidae" "Paralichthys"
Seq_7 "Eukaryota" "Chordata" "Actinopteri" "Clupeiformes" "Clupeidae" "Alosa"
Seq_9 "Eukaryota" "Chordata" "Actinopteri" "Gobiiformes" "Gobiidae" "Gobiosoma"
Seq_10 "Eukaryota" "Chordata" "Actinopteri" "Pleuronectiformes" "Scophthalmidae" "Scophthalmus"
Seq_11 "Eukaryota" "Chordata" "Actinopteri" "Perciformes" "Serranidae" "Centropristis"
Seq_12 "Eukaryota" "Chordata" "Actinopteri" "Spariformes" "Sparidae" "Stenotomus"
Seq_15 "Eukaryota" "Chordata" "Actinopteri" NA "Sciaenidae" "Leiostomus"
Seq_16 "Eukaryota" "Chordata" "Actinopteri" NA "Sciaenidae" "Menticirrhus"
Seq_17 "Eukaryota" "Chordata" "Actinopteri" "Labriformes" "Labridae" "Tautoga"
Seq_19 "Eukaryota" "Chordata" "Actinopteri" "Perciformes" "Cottidae" "Myoxocephalus"
Seq_20 "Eukaryota" "Chordata" "Actinopteri" "Pleuronectiformes" "Pleuronectidae" "Pseudopleuronectes"
Seq_21 "Eukaryota" "Chordata" "Actinopteri" NA "Moronidae" "Morone"
Seq_22 "Eukaryota" "Chordata" "Actinopteri" "Syngnathiformes" "Syngnathidae" "Syngnathus"
Seq_30 "Eukaryota" "Chordata" "Actinopteri" "Pleuronectiformes" "Paralichthyidae" "Etropus"
Seq_33 "Eukaryota" "Chordata" "Actinopteri" NA "Sciaenidae" "Cynoscion"
Seq_34 "Eukaryota" "Chordata" "Actinopteri" "Labriformes" "Labridae" "Tautogolabrus"
Seq_36 "Eukaryota" "Chordata" "Actinopteri" "Anguilliformes" "Anguillidae" "Anguilla"
Seq_38 "Eukaryota" "Chordata" "Actinopteri" "Scombriformes" "Scombridae" "Thunnus"
Seq_40 "Eukaryota" "Chordata" "Actinopteri" "Perciformes" "Gasterosteidae" "Apeltes"
Seq_44 "Eukaryota" "Chordata" "Actinopteri" "Cyprinodontiformes" "Fundulidae" "Fundulus"
Seq_50 "Eukaryota" "Chordata" "Actinopteri" "Atheriniformes" "Atherinopsidae" "Membras"
Seq_52 "Eukaryota" "Chordata" "Actinopteri" "Gadiformes" "Phycidae" "Urophycis"
Seq_54 "Eukaryota" "Chordata" "Actinopteri" "Scombriformes" "Scombridae" "Scomber"
Seq_57 "Eukaryota" "Chordata" "Actinopteri" "Perciformes" "Triglidae" "Prionotus"
Seq_67 "Eukaryota" "Chordata" "Actinopteri" "Scombriformes" "Scombridae" "Thunnus"
Seq_82 "Eukaryota" "Chordata" "Actinopteri" NA "Sciaenidae" "Bairdiella"
Seq_84 "Eukaryota" "Chordata" "Actinopteri" "Gadiformes" "Gadidae" "Microgadus"
Seq_102 "Eukaryota" "Chordata" "Actinopteri" "Clupeiformes" "Engraulidae" "Engraulis"
Seq_103 "Eukaryota" "Chordata" "Actinopteri" "Perciformes" "Cottidae" "Myoxocephalus"
Seq_115 "Eukaryota" "Chordata" "Actinopteri" "Cyprinodontiformes" "Fundulidae" "Fundulus"
Seq_139 "Eukaryota" "Chordata" "Actinopteri" "Batrachoidiformes" "Batrachoididae" "Opsanus"
Seq_141 "Eukaryota" "Chordata" "Actinopteri" "Scombriformes" "Scombridae" "Katsuwonus"
Seq_181 "Eukaryota" "Chordata" "Actinopteri" "Tetraodontiformes" "Tetraodontidae" "Sphoeroides"
Seq_231 "Eukaryota" "Chordata" "Actinopteri" "Gadiformes" "Merlucciidae" "Merluccius"
Seq_359 "Eukaryota" "Chordata" "Actinopteri" "Perciformes" "Triglidae" "Prionotus"
Seq_362 "Eukaryota" "Chordata" "Chondrichthyes" "Myliobatiformes" "Myliobatidae" "Rhinoptera"
Seq_372 "Eukaryota" "Chordata" "Chondrichthyes" "Carcharhiniformes" "Triakidae" "Mustelus"
species
Seq_1 "Menidia menidia"
Seq_2 "Brevoortia tyrannus"
Seq_3 "Anchoa mitchilli"
Seq_4 "Pomatomus saltatrix"
Seq_5 "Lutjanus griseus"
Seq_6 "Paralichthys dentatus"
Seq_7 "Alosa mediocris"
Seq_9 "Gobiosoma ginsburgi"
Seq_10 "Scophthalmus aquosus"
Seq_11 "Centropristis striata"
Seq_12 "Stenotomus chrysops"
Seq_15 "Leiostomus xanthurus"
Seq_16 "Menticirrhus saxatilis"
Seq_17 "Tautoga onitis"
Seq_19 "Myoxocephalus aenaeus"
Seq_20 "Pseudopleuronectes americanus"
Seq_21 "Morone saxatilis"
Seq_22 "Syngnathus fuscus"
Seq_30 "Etropus microstomus"
Seq_33 "Cynoscion regalis"
Seq_34 "Tautogolabrus adspersus"
Seq_36 "Anguilla rostrata"
Seq_38 "Thunnus obesus"
Seq_40 "Apeltes quadracus"
Seq_44 "Fundulus majalis"
Seq_50 "Membras martinica"
Seq_52 "Urophycis floridana"
Seq_54 "Scomber japonicus"
Seq_57 "Prionotus carolinus"
Seq_67 "Thunnus thynnus"
Seq_82 "Bairdiella chrysoura"
Seq_84 "Microgadus tomcod"
Seq_102 "Engraulis mordax"
Seq_103 "Myoxocephalus quadricornis"
Seq_115 "Fundulus heteroclitus"
Seq_139 "Opsanus tau"
Seq_141 "Katsuwonus pelamis"
Seq_181 "Sphoeroides maculatus"
Seq_231 "Merluccius bilinearis"
Seq_359 "Prionotus evolans"
Seq_362 "Rhinoptera bonasus"
Seq_372 "Mustelus canis"
NOTES for Sara
I am actually getting 42 unique species- which ones am I missing that should be removed?
Based on my previous scripts with Cariaco Eukaryotic data
# convert ps object to dataframe using phyloseq's psmelt
species_df <- psmelt(speciesGlommed_RA)
# replace zeroes in the table with NA
species_df[species_df == 0] <- NA
# and remove rows with NAs in abundance (this is so they don't appear as small dots in plot)
species_df <- filter(species_df, !is.na(Abundance))
speciesbubbleplot <- ggplot(species_df, aes(x = Station, y = fct_rev(species), color = Station)) + # the fancy stuff around y (species) helps to present it in reverse order in the plot (from top to btm alphabetically)
geom_point(aes(size = Abundance, fill = Station), color = "black", pch = 21)+
scale_size(range = c(1,15)) +
scale_size_area(breaks = c(0,.25,.5,.75,1), max_size = 6)+
xlab("")+
ylab("")+
labs(size="Relative Abundance")+
theme_bw() +
scale_fill_brewer(palette="Paired") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank()) +
facet_grid(Datecode~Bayside, scales = "free", space = "free", drop= TRUE)
Scale for 'size' is already present. Adding another scale for 'size', which will replace the existing scale.
speciesbubbleplot
I like the layout above so I will export
ggsave(filename = "Figures/speciesbubbleplot.eps", plot = speciesbubbleplot, units = c("in"), width = 7, height = 12, dpi = 300, )
NOTES next Sara wants the average rel abundances for each species across all dates
PCA is essentially a type of PCoA using the Euclidean distance matrix as input. When combined with a log-ratio transformation of the count table, this is deemed appropriate for compositional datasets. It is also recommended as a first step in exploratory analyses.
First do a CLR, centered log ratio transformation of the absolute abundance data (after filtering), as suggested by Gloor et al. 2017
# Estimate covariance matrix for CLR-transformed ASV table
clr_asv_table_ps <- data.frame(compositions::clr(otu_table(ps)))
Generate the PCA and visualize axes
# Generate a Principle Component Analysis (PCA) and evaluated based on the eigen decomposition from sample covariance matrix.
lograt_pca <- prcomp(clr_asv_table_ps)
# NOTE- this is equivalent to first making a Euclidean distance matrix using the CLR data table and then running a PCoA. A Euclidean distance matrix of a log-transformed data table = an Aitchison distance matrix. So this is equivalent to the compositional methods listed in Gloor et al.
# Visual representation with a screeplot
lograt_variances <- as.data.frame(lograt_pca$sdev^2/sum(lograt_pca$sdev^2)) %>% #Extract axes
# Format to plot
select(PercVar = 'lograt_pca$sdev^2/sum(lograt_pca$sdev^2)') %>%
rownames_to_column(var = "PCaxis") %>%
data.frame
head(lograt_variances)
# Plot screeplot
ggplot(lograt_variances, aes(x = as.numeric(PCaxis), y = PercVar)) +
geom_bar(stat = "identity", fill = "grey", color = "black") +
theme_minimal() +
theme(axis.title = element_text(color = "black", face = "bold", size = 10),
axis.text.y = element_text(color = "black", face = "bold"),
axis.text.x = element_blank()) +
labs(x = "PC axis", y = "% Variance", title = "Log-Ratio PCA Screeplot, CLR Tranformation")
A faithful representation is to plot this ordination in 2D because the 3rd and 4th axes explain very low % of variances, while the 1st and 2nd explain a decently large proportion of variance (15.7 + 10.5 = 26.2%)
Visualize the PCA-
lograt_pca$x #View PC values
PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
T1S10 -2.6238954 -0.3252536 1.27005712 -1.666211452 0.798649968 -0.57848189 0.247580998 -0.558711581 1.54051341 -1.00257763
T1S11 0.2705609 -4.0749320 -6.66931720 2.284305153 2.672326913 -0.41018048 -0.526654868 -0.462324083 -0.41629925 -1.22992687
T1S1 -2.7027302 2.5315244 0.42770674 -2.891886434 -0.089199476 -2.02741275 -1.526307265 4.507918543 -2.80161601 0.99281647
T1S3 -2.1044928 1.0756380 1.49045892 -2.219957763 -0.061377153 -0.38494915 0.320469272 -1.012257490 1.42749170 -0.65295542
T1S5 -1.7135917 1.1790076 2.77474051 -1.340062932 -1.364232051 0.24102503 -0.009868407 -0.958600452 1.66240096 1.90249139
T1S6 -1.6263823 -0.7382611 4.92533402 -2.100496346 -0.992855481 1.31005659 2.080679436 -0.291113329 -0.53023688 -0.48181295
T1S9 -1.8890736 -7.1882980 0.74791280 0.478128606 -1.223416026 -1.18521711 2.834417355 1.011107634 -1.91733560 -1.95236999
T2S10 -5.4723849 -3.7102032 -2.02733221 0.358739428 0.095792276 -4.19640555 1.696323226 0.708131115 1.11063176 -1.77725044
T2S11 4.0600493 -1.6653347 0.95792661 -1.382521282 -1.935788288 0.14614516 -0.219686405 0.066454230 0.28910892 -0.35381697
T2S1 -0.5256139 -1.6043596 -2.82111238 0.471370176 5.396921155 4.13619608 0.992107496 -1.578763421 1.80377493 2.18025904
T2S2 -1.5675813 3.4039386 -2.83994788 -1.328904319 1.293026359 3.90586508 4.611348013 -0.711643119 -2.75478280 1.19354776
T2S3 -2.4949331 1.8527466 0.55242812 -2.302859249 0.001024354 -1.64228999 -0.599524728 -0.235068104 0.87739077 -0.06120956
T2S4 -5.4192303 -0.5712951 -2.80407160 0.328311261 -0.473384127 3.55256896 -0.864999520 3.421745399 -2.94070397 -1.64572626
T2S5 -2.2024924 2.1002160 2.34821607 2.710381042 3.191568957 -0.90891012 -3.088217336 3.454337657 0.47494384 2.76967314
T2S6 -1.8637965 -5.0301644 2.83044061 0.651949689 2.536239781 2.28530167 -1.412309984 2.720676493 2.49919366 1.43936846
T2S9 -8.1183606 -4.0795361 -5.89618438 3.088417424 -6.620641957 0.18089576 0.254590155 -0.882108849 1.90441468 4.35599519
T3S10 4.6214005 -1.6575522 0.58946754 -1.320385118 -1.661933332 0.18164541 -0.035201214 0.001341068 0.26754996 -0.43305903
T3S1 -2.3913351 2.2896417 -0.20957141 -1.068738089 -0.640139573 0.63357839 -1.285553135 -1.195502284 -0.53495156 -0.32737072
T3S2 -2.2678012 3.3894600 0.27244730 2.244919808 0.035628741 0.18047505 -3.537238882 -0.264238101 -0.35901011 0.87134641
T3S3 -0.6141891 1.8274067 -1.15847880 -1.437042040 0.895015921 0.69099724 2.030728395 -1.239263664 -0.42837700 0.42179666
T3S4 -3.4730868 0.5900791 0.55663668 -1.944078098 1.248042434 -0.03954649 0.588369559 -0.373556187 3.33127230 -1.17587526
PC11 PC12 PC13 PC14 PC15 PC16 PC17 PC18 PC19 PC20
T1S10 0.22597717 -0.49591596 0.62088813 -0.11135149 -0.90317596 0.052175005 0.51636597 0.11646915 -1.13617561 -0.939740199
T1S11 -0.33875912 -2.47937788 -1.92292360 -1.61727672 1.42629776 2.633399673 -1.77511843 -0.64393857 -0.16379258 -1.170843742
T1S1 0.42101859 0.65442210 0.05542544 1.05706121 1.01377542 1.748809810 0.60333600 -2.11186111 -3.09626558 1.781651311
T1S3 0.21932314 -0.65001251 0.19899498 -0.42992913 -0.17187489 -0.321107836 -0.29559955 -0.62907644 -0.98165385 -0.420964761
T1S5 0.90236276 2.37041721 -0.90409728 -5.95846889 -0.99039907 2.323748019 1.87406642 -0.52235594 0.47615803 -0.538955197
T1S6 0.32492687 -0.58250914 -3.36962676 1.03850812 2.92510401 0.527322477 0.74283517 1.98271154 -0.90769777 0.208257905
T1S9 -1.82706149 4.63700830 -1.41096275 1.65684273 -0.73128430 -0.439014716 0.22456025 -2.83124825 1.54026880 -1.804244119
T2S10 -0.63570055 -0.04821840 0.03027849 -2.11585154 1.52888813 -1.711272217 -2.05352206 -0.32604955 -0.43090515 1.805775529
T2S11 0.10750821 -0.16232384 0.22061834 0.06591853 -0.35834169 -0.180707353 -0.42978954 0.20962610 -0.45129627 -0.074319773
T2S1 -0.84361171 2.61168296 -1.64637670 0.47643508 -1.53112443 -0.783075685 1.14544456 0.43919900 -1.03784578 1.943550414
T2S2 1.90947454 -0.88042217 2.50315503 -0.71782729 2.65945940 -0.679764282 0.20848879 -1.61984898 0.76279443 -1.652676887
T2S3 -0.28832020 -0.58979394 0.08413418 -0.03476737 0.53372457 -0.528477113 -0.10571160 -0.13196257 -0.50579456 -0.513744834
T2S4 -1.73341212 -0.05444488 2.67180984 -1.40242616 -0.29205849 0.287767338 0.90589382 2.03117571 -0.08582432 0.517562229
T2S5 2.20888695 1.28719261 0.16320560 0.97886685 -0.43996377 -0.185431085 0.30838611 -0.78176830 0.76544762 -0.767696679
T2S6 -2.05629428 -1.67779276 0.06090643 -0.37825238 1.14814511 0.229027684 -1.04539827 -0.09822133 0.79124945 -0.201171927
T2S9 2.52542213 -1.21424602 -0.33127400 1.97232751 -0.63646116 0.587380408 0.56256483 0.20782020 0.12040843 0.385932778
T3S10 0.17608809 -0.09239084 0.27696903 0.01189034 -0.04774737 -0.240905695 -0.07265447 0.14762881 -0.33514522 0.142036307
T3S1 -1.77915713 -1.57684401 -2.40731062 -0.16246954 -1.13827510 -1.572223128 0.39958868 -0.62003482 -1.23180650 -0.908141905
T3S2 -1.80093781 1.40677342 -2.00653035 1.31217445 2.49078997 -0.331605766 -0.44636779 1.88665314 0.62390115 -1.979127915
T3S3 -0.19515396 0.64121705 -1.42175886 0.94633601 1.18314322 0.356666826 0.35724747 -0.39248176 -1.01111108 0.887210503
T3S4 -0.03337711 -1.22064722 1.41331871 0.55048179 -0.99424565 -0.090116768 -0.12121385 -1.12841112 0.14705172 -0.272828981
PC21 PC22 PC23 PC24 PC25 PC26 PC27 PC28 PC29
T1S10 0.14553176 0.444887701 0.030555073 -0.10425935 -0.4649432055 0.355881811 0.17992733 0.7760907427 -0.627690690
T1S11 1.38300176 0.979099117 0.211703218 0.15116950 -0.8685434581 1.049405232 -1.10236127 0.0002552213 -0.767570144
T1S1 0.34009906 -0.271639014 -0.903345414 1.30069377 -0.5473422995 -0.048022311 -0.10948686 -0.9587529174 -0.015703705
T1S3 0.86994683 0.246913911 -0.158519509 0.12505746 -1.1016182738 0.059396108 0.01323105 1.2624136272 0.275632635
T1S5 -0.28004625 -0.600135978 -0.225638160 -0.22408151 -0.0009001045 -0.794429743 -0.71040916 -0.7007108564 -0.013510170
T1S6 0.30251073 1.892988844 1.479249040 -2.05668429 0.6892908892 -0.364245835 -0.08183022 -0.8163892316 0.653306565
T1S9 0.60401701 -0.026577052 -0.266982730 -0.05322975 0.3969389259 0.099425676 -0.35147954 0.1775020985 0.456170749
T2S10 -2.33247107 -1.423626311 0.484900870 -1.55699001 -0.7921676214 -0.407380274 0.94742076 0.0174250216 -0.392227697
T2S11 -0.37304090 -0.044363953 -0.212358865 0.03656853 -0.2055835293 0.405011932 -0.12423554 -0.1260813064 0.255872703
T2S1 0.51211253 -0.869077708 0.215385184 -0.23843593 -1.4081013750 1.359543178 0.89307294 -0.3493111760 0.367630027
T2S2 -1.23341832 0.332054911 -0.742253664 0.06444234 0.3595826030 0.578803687 1.26134666 -0.7590524011 -0.330851226
T2S3 1.15014708 0.327327365 -0.825669816 -0.49833503 -0.5864129208 0.557311711 0.31739992 1.1282591094 0.969019040
T2S4 0.65239653 -0.397745489 -0.437515949 -1.40322970 0.5273097762 0.006654982 -1.03730521 1.4003376711 0.686279185
T2S5 -1.53158603 1.099204218 0.981168725 -1.53870366 0.6821822418 0.606535859 -0.09598065 1.0952280391 -0.572618674
T2S6 1.14916912 -0.442956604 -0.031343541 1.53126698 1.1626853043 -2.294630143 1.30790014 -0.4321071978 0.039003479
T2S9 -0.06879790 0.314833969 -0.487910272 0.13571852 -0.0935417744 -0.111592202 0.13245163 0.0115432579 0.137289451
T3S10 -0.31137995 0.003053146 -0.131364000 0.14127153 -0.1147585816 0.378698308 -0.11600723 0.0008098405 -0.054375007
T3S1 -0.18259033 -0.941533760 -0.608857056 -0.32371094 1.3136321089 -0.314823572 1.28871003 0.3636956219 0.116792357
T3S2 -1.90095709 -1.656661333 -0.169026823 1.05756571 -0.7505352043 0.531616513 -1.01059260 -0.4445991478 -0.217363962
T3S3 0.23648113 0.150139986 0.744726257 1.08530156 0.9973389162 -0.297207577 -0.70397225 1.0798450685 -0.540203211
T3S4 -1.40747880 0.950362687 2.174998385 1.84091769 0.5522902416 0.509635655 -0.78582948 0.1257947674 1.529647938
PC30 PC31 PC32 PC33 PC34 PC35 PC36 PC37 PC38
T1S10 -0.61090108 -0.29334392 -0.435556282 0.436122013 0.620570063 0.0752049091 0.364291442 0.222722389 0.853366974
T1S11 0.38683053 -0.47170713 -0.507376193 -0.200978899 -0.466127983 -0.0887001046 0.407399148 0.137577782 -0.114942713
T1S1 0.39068097 -0.29859127 0.422021427 0.016638837 0.130295349 -0.3423059226 0.099331214 0.004371715 -0.053692539
T1S3 -0.89085871 0.30772799 0.035796223 0.384888434 1.245611088 0.8423535329 -0.056811384 0.677715223 0.876394557
T1S5 -0.10448823 -0.34639788 -0.108416705 0.280931290 -0.214231135 0.2227704635 0.234041079 0.027794099 -0.168256122
T1S6 0.17543435 -0.16388918 0.408543638 -0.023459302 -0.021177178 -0.1444244416 0.084570850 0.490213102 -0.165936538
T1S9 -0.06692976 0.30760983 -0.009473703 0.040284195 -0.050039551 0.0235304169 0.049103401 0.134936707 -0.032580927
T2S10 0.45881864 -0.27229696 0.378141816 0.079696113 0.025621644 0.2220459377 -0.229916233 -0.116912242 0.023235021
T2S11 -0.07367150 0.04060642 0.327821556 -0.335943407 0.177567813 -0.2684558422 0.561488607 0.476894504 0.234251029
T2S1 0.22719978 0.38827525 0.282206264 -0.438390883 0.036965892 -0.1441772303 0.141507739 0.073251632 -0.131871600
T2S2 -0.33864908 0.16809097 0.064604371 0.086877729 0.002553673 0.1239108829 0.412985888 0.139172474 0.005935413
T2S3 -2.26654726 0.23001824 0.880133582 -0.079399503 -1.611846710 -0.0529833199 -0.104266043 -0.665332121 -0.874030961
T2S4 0.77112458 0.73081680 -0.074041702 0.112370407 0.295223675 0.4241879490 -0.004228575 0.091313711 -0.312300967
T2S5 -0.17121290 -0.74735265 -0.331455620 -0.613891167 -0.362680848 -0.2437333112 -0.107411964 0.564034865 0.585280991
T2S6 -0.45836363 0.59828788 0.257913281 -0.377057458 -0.112145871 0.2661286081 0.001694725 0.444611086 0.313382173
T2S9 0.04952467 0.28054101 -0.020026648 -0.062809756 0.036918082 0.0670041757 -0.061020750 -0.038346173 0.043463127
T3S10 -0.04433020 0.11256207 0.276748498 0.012290749 0.105555296 -0.3108955253 0.225398056 0.383416618 0.104263911
T3S1 0.78849208 -0.09561134 -1.600916711 1.062468681 -0.458306008 -1.1614572080 0.283387190 -0.270171409 0.076478396
T3S2 -0.40605694 0.55137855 0.461720824 0.721368991 0.554321865 0.4202858684 0.512499024 -0.227527215 -0.125192170
T3S3 0.02013393 -0.20179302 -0.865030680 0.385797223 -0.040604332 1.1497492465 -1.231478661 -0.824389167 0.413980530
T3S4 0.97160611 -0.31764671 0.039652545 0.203106279 0.405656989 -0.0564320934 0.292736573 0.165312522 -1.225094311
PC39 PC40 PC41 PC42 PC43 PC44 PC45 PC46 PC47
T1S10 -0.38601496 -0.522398117 -0.2601828761 1.7425428517 -0.150037940 -0.022197976 -2.831425e-03 -0.0100519518 -2.038300e-16
T1S11 -0.22022127 -0.233355749 -0.0045061584 -0.2088644668 -0.045195931 0.051624214 6.578873e-03 0.0074415753 2.957704e-16
T1S1 -0.21897172 -0.033400842 0.1854388730 0.1027272364 0.079799837 -0.025304143 -5.997196e-03 0.0027479728 -9.393528e-16
T1S3 -0.75497648 0.603862482 0.6040479086 -1.0356859876 0.297725332 -0.016166585 -1.318602e-02 0.0043860796 -1.147520e-15
T1S5 0.26017909 0.053160901 -0.0241782401 -0.0575928258 0.020748744 -0.009410090 9.226674e-05 -0.0042139111 2.766017e-15
T1S6 -0.25843482 0.293010354 0.0953476162 0.1450311162 0.083660353 0.004693246 -1.979926e-02 0.0067988303 3.126839e-15
T1S9 -0.19310548 0.050539539 -0.0263737418 0.0752063777 0.060796962 0.005479383 -4.041996e-03 0.0066798004 5.663872e-16
T2S10 0.30388992 0.100877709 -0.0191138984 0.0412173967 0.026576122 -0.032226350 -3.521885e-03 -0.0070238059 -3.287301e-16
T2S11 0.28580250 -0.321592098 -0.4781464260 -0.4109621304 -1.159042956 -1.198398284 -2.007390e-01 0.0455439782 -1.189153e-15
T2S1 -0.11248568 -0.073344813 -0.0955302997 -0.0376190978 -0.004465990 0.023019713 8.557707e-03 0.0032389114 -1.001803e-15
T2S2 0.01410303 -0.018399653 -0.0723573078 -0.0394075311 0.055382093 -0.053886132 2.126331e-02 -0.0041445615 1.569925e-16
T2S3 0.47491853 -0.324154290 0.1113297600 0.0257086927 0.021247549 -0.026144440 -1.289991e-03 -0.0004957201 -1.286297e-15
T2S4 0.02175802 -0.371533083 -0.0654592950 -0.0335770823 -0.032892383 0.006709334 4.995619e-03 -0.0041659031 3.790371e-16
T2S5 0.28095908 -0.242536910 0.0130502993 -0.3851088178 -0.060000171 0.091145738 9.202016e-03 0.0107147824 7.953707e-16
T2S6 0.17245826 -0.168809317 -0.0903469935 0.0539771988 -0.089777526 0.067356765 4.954736e-03 0.0027736662 2.391316e-15
T2S9 -0.04763513 -0.036926578 0.0312385396 0.0438364840 0.023883432 0.018597761 5.213981e-03 0.0019178143 -2.351418e-15
T3S10 0.34199375 -0.272932684 -0.2206710864 -0.2167177793 -0.243551170 0.556300282 2.557553e-01 -0.8175569629 -5.924081e-16
T3S1 -0.17851876 -0.315320817 0.4560668700 -0.3531386713 -0.169198752 0.016921196 -1.278480e-02 -0.0068107404 -9.254750e-16
T3S2 0.12952809 -0.070809741 0.1503751764 0.0970471976 0.009140306 -0.010012191 -6.877311e-03 -0.0104988319 -3.981190e-16
T3S3 1.09009743 0.242828505 -0.7044769588 -0.1231401097 -0.397445764 -0.026921698 2.386564e-02 -0.0229946682 -2.870967e-16
T3S4 0.33863712 -0.189836596 0.1949792205 0.0225704339 0.229417048 -0.063472056 -6.409916e-03 0.0074190104 1.847481e-16
[ reached getOption("max.print") -- omitted 26 rows ]
Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls
4: In readChar(file, size, TRUE) : truncating string with embedded nuls
5: In readChar(file, size, TRUE) : truncating string with embedded nuls
6: In readChar(file, size, TRUE) : truncating string with embedded nuls
7: In readChar(file, size, TRUE) : truncating string with embedded nuls
8: In readChar(file, size, TRUE) : truncating string with embedded nuls
9: In readChar(file, size, TRUE) : truncating string with embedded nuls
10: In readChar(file, size, TRUE) : truncating string with embedded nuls
pca_lograt_frame <- data.frame(lograt_pca$x) %>%
rownames_to_column(var = "SampleID")
# Merge metadata into the pca data table
pca_lograt_frame <- left_join(pca_lograt_frame, metadata, by = "SampleID")
head(pca_lograt_frame)
# Plot PCA with Station
pca_lograt_station <- ggplot(pca_lograt_frame, aes(x = PC1, y = PC2, color = Station)) +
geom_point(aes(shape = Bayside), size = 4) +
ylab(paste0('PC2 ', round(lograt_variances[2,2]*100,2),'%')) + #Extract y axis value from variance
xlab(paste0('PC1 ', round(lograt_variances[1,2]*100,2),'%')) + #Extract x axis value from variance
scale_color_brewer(palette="Paired") +
ggtitle('CLR-Euclidean PCA') +
coord_fixed(ratio = 1) +
theme_bw()
pca_lograt_station
ggsave("figures/pca_clr.eps",pca_lograt_station, width = 7, height = 5, units = c("in"))
The CLR-Euclidean PCA reveals there is some separation according to East vs West. The PCA only explains ~26% of the variance so keep going with different ordinations to see if we can get a better representation
The more traditional approach to ordinations is to do a PCoA on a distance matrix such as Bray-Curtis, Jaccard, or Unifrac. While these are not considered compositional approaches, when combined with pre-treatment (transformations) they become more appropriate. One such transformation that I will use here is the Hellinger transformation.
The different distance matrices also tell you a few different things about the dataset so I will run through this to try to see if I can tease those out.
Before calculating any distance matrix, do a transformation of the filtered count table. Hellinger transformation is the square root of the relative abundance, so calculate it based on the ps_ra object:
ps_hellinger
phyloseq-class experiment-level object
otu_table() OTU Table: [ 368 taxa and 47 samples ]
sample_data() Sample Data: [ 47 samples by 8 sample variables ]
tax_table() Taxonomy Table: [ 368 taxa by 7 taxonomic ranks ]
First, Jaccard, which builds the distance matrix based on presence/absence between samples. It does not take into account relative abundance of the taxa. Therefore this functions well for determining differences driven by rare taxa, which are weighed the same as abundant taxa.
jac_dmat<-vegdist(otu_table(ps_hellinger),method="jaccard") # Jaccard dist metric
pcoa_jac<-pcoa(jac_dmat) # perform PCoA
# Extract variances from pcoa, from jaccard calculated dist. metric
jac_variances <- data.frame(pcoa_jac$values$Relative_eig) %>%
select(PercVar = 'pcoa_jac.values.Relative_eig') %>%
rownames_to_column(var = "PCaxis") %>%
data.frame
head(jac_variances)
# Make a screeplot
ggplot(jac_variances, aes(x = as.numeric(PCaxis), y = PercVar)) +
geom_bar(stat = "identity", fill = "grey", color = "black") +
theme_minimal() +
theme(axis.title = element_text(color = "black", face = "bold", size = 10),
axis.text.y = element_text(color = "black", face = "bold"),
axis.text.x = element_blank()) +
labs(x = "PC axis", y = "% Variance", title = "Jaccard PCoA Screeplot")
The first two axes (19.0 + 9.6 = 28.6) are OK. But I am going to experiment and plot the first 3 axes since the 2nd and 3rd explain a similar amount of variance, 19.6 and 8.4% respectively
Plot in 3D with Plotly
# Extract variances from the jaccard pcoa
pcoa_jac_df <- data.frame(pcoa_jac$vectors) %>%
rownames_to_column(var = "SampleID")
# Merge metadata into the pcoa data table
pcoa_jac_df <- left_join(pcoa_jac_df, metadata, by = "SampleID")
head(pcoa_jac_df)
# Select eigenvalues from dataframe, round to 4 places and multiply by 100 for plotting. These will be the axes for the 3-D plot
eigenvalues<-round(jac_variances[,2], digits = 4)*100
# Plotly - 3-D
pcoa_jaccard <- plot_ly(pcoa_jac_df, type='scatter3d', mode='markers',
x=~Axis.2,y=~Axis.3,z=~Axis.1,colors=~brewer.pal(11,'Paired'),
color=~Station, symbols = c('circle','diamond'), symbol=~Bayside)%>%
layout(font=list(size=12),
title='PCoA Jaccard Distance',
scene=list(xaxis=list(title=paste0('Co 2 ',eigenvalues[2],'%'),
showticklabels=FALSE,zerolinecolor='black'),
yaxis=list(title=paste0('Co 3 ',eigenvalues[3],'%'),
showticklabels=FALSE,zerolinecolor='black'),
zaxis=list(title=paste0('Co 1 ',eigenvalues[1],'%'),
showticklabels=FALSE,zerolinecolor='black')))
pcoa_jaccard
withr::with_dir('Figures', htmlwidgets::saveWidget(as_widget(pcoa_jaccard), file="pcoa_jaccard.html"))
The Jaccard-PCoA shows separation along axis 2 in East vs West differences.
Next, try a Bray-Curtis distance matrix with PCoA, which builds the distance matrix based on presence/absence between samples and relative abundance differences. This ordination will represent well the differences in samples that are driven by taxa with high relative abundances.
bray_dmat<-vegdist(otu_table(ps_hellinger),method="bray") # Bray-Curtis dist metric
pcoa_bray<-pcoa(bray_dmat) # perform PCoA
# Extract variances from pcoa, from jaccard calculated dist. metric
bray_variances <- data.frame(pcoa_bray$values$Relative_eig) %>%
select(PercVar = 'pcoa_bray.values.Relative_eig') %>%
rownames_to_column(var = "PCaxis") %>%
data.frame
head(bray_variances)
# Make a screeplot
ggplot(bray_variances, aes(x = as.numeric(PCaxis), y = PercVar)) +
geom_bar(stat = "identity", fill = "grey", color = "black") +
theme_minimal() +
theme(axis.title = element_text(color = "black", face = "bold", size = 10),
axis.text.y = element_text(color = "black", face = "bold"),
axis.text.x = element_blank()) +
labs(x = "PC axis", y = "% Variance", title = "Bray-Curtis PCoA Screeplot")
The first two axes (27.7 + 14.3 = 42%) are pretty good again but I am still going to experiment in the plot with the 3rd axis since it is similar to the second (12.2% variance)
Plot in 3D with Plotly
# Extract variances from the jaccard pcoa
pcoa_bray_df <- data.frame(pcoa_bray$vectors) %>%
rownames_to_column(var = "SampleID")
# Merge metadata into the pcoa data table
pcoa_bray_df <- left_join(pcoa_bray_df, metadata, by = "SampleID")
head(pcoa_bray_df)
# Select eigenvalues from dataframe, round to 4 places and multiply by 100 for plotting. These will be the axes for the 3-D plot
eigenvalues<-round(bray_variances[,2], digits = 4)*100
# Plotly - 3-D
pcoa_bray <- plot_ly(pcoa_bray_df, type='scatter3d', mode='markers',
x=~Axis.2,y=~Axis.3,z=~Axis.1,colors=~brewer.pal(11,'Paired'),
color=~Station, symbols = c('circle','diamond'), symbol=~Bayside)%>%
layout(font=list(size=12),
title='PCoA Bray-Curtis Distance',
scene=list(xaxis=list(title=paste0('Co 2 ',eigenvalues[2],'%'),
showticklabels=FALSE,zerolinecolor='black'),
yaxis=list(title=paste0('Co 3 ',eigenvalues[3],'%'),
showticklabels=FALSE,zerolinecolor='black'),
zaxis=list(title=paste0('Co 1 ',eigenvalues[1],'%'),
showticklabels=FALSE,zerolinecolor='black')))
pcoa_bray
withr::with_dir('Figures', htmlwidgets::saveWidget(as_widget(pcoa_bray), file="pcoa_bray.html"))
These results are similar to Jaccard: the second axis seems driven by differences in East vs West. But there are clearly other things going on here with axies 1 and 3. I think this is a good representation of the data: together the 3 axes explain 54.13% of the variance.
Lastly, try a non-metric dimensional scaling ordination. PCA/PCoA are metric and attempt to rotate axes to fit the distance matrix distribution. An NMDS represents the data in 2-axes, by constraining the distribution of the points. Similar to above, this can be combined with different pre-treatment of the data.
First try the compositional approach, an NMDS on CLR-tranformed data using the Euclidean distances (aka Aitchison distance)
euc_dmat<-dist(clr_asv_table_ps, method = "euclidean") # Build the Aitchison distance matrix
euc_nmds <- metaMDS(euc_dmat, k=2, autotransform=FALSE) # Run the ordination
Run 0 stress 0.2105436
Run 1 stress 0.2128996
Run 2 stress 0.2122615
Run 3 stress 0.2111085
Run 4 stress 0.2130327
Run 5 stress 0.223828
Run 6 stress 0.2127398
Run 7 stress 0.2171268
Run 8 stress 0.2233421
Run 9 stress 0.2390473
Run 10 stress 0.2232293
Run 11 stress 0.2146083
Run 12 stress 0.2343976
Run 13 stress 0.2122711
Run 14 stress 0.2175982
Run 15 stress 0.2273345
Run 16 stress 0.227759
Run 17 stress 0.2231565
Run 18 stress 0.2286921
Run 19 stress 0.2426823
Run 20 stress 0.2110142
... Procrustes: rmse 0.01362479 max resid 0.06588342
*** No convergence -- monoMDS stopping criteria:
1: no. of iterations >= maxit
19: stress ratio > sratmax
euc_nmds$stress #Check the stress. Less than 0.1 is good. Less than 0.05 is better. This will be different each time, since it is iteratively finding a unique solution each time (although the should look similar)
[1] 0.2105436
# Extract points from nmds and merge into data frame with metadata
euc_nmds_df <- data.frame(euc_nmds$points) %>%
rownames_to_column(var = "SampleID")
# Merge metadata into the pcoa data table
euc_nmds_df <- left_join(euc_nmds_df, metadata, by = "SampleID")
head(euc_nmds_df)
## Plotting euclidean distance NMDS
nmds_aitch <- ggplot(euc_nmds_df,aes(x = MDS1, y = MDS2, color = Station, shape = Bayside)) +
geom_point(size = 4) +
scale_color_brewer(palette="Paired") +
theme_bw() +
labs(x = "NMDS 1", y = "NMDS 2", title = paste0('Aitchison Distance NMDS, Stress = ', round(euc_nmds$stress,2))) +
coord_fixed(ratio = 1)
nmds_aitch
ggsave("figures/nmds_aitch.eps",nmds_aitch, width = 7, height = 5, units = c("in"))
The above has a relatively high stress (>0.2) so should be interpreted with caution. But it does show some separation East vs West along NMDS 1.
Next try a Jaccard NMDS, which will represent differences in presence/absence among samples, emphasizing both abundant and rare taxa the same
jac_nmds <- metaMDS(jac_dmat, k=2, autotransform=FALSE) # Run the ordination. Distance matrix was already calculated above
Run 0 stress 0.1627003
Run 1 stress 0.1660244
Run 2 stress 0.1496856
... New best solution
... Procrustes: rmse 0.09100415 max resid 0.3183697
Run 3 stress 0.1495161
... New best solution
... Procrustes: rmse 0.05088368 max resid 0.3258651
Run 4 stress 0.1578622
Run 5 stress 0.1664404
Run 6 stress 0.1572462
Run 7 stress 0.1573408
Run 8 stress 0.1495163
... Procrustes: rmse 0.0009228226 max resid 0.004331749
... Similar to previous best
Run 9 stress 0.1496491
... Procrustes: rmse 0.01267127 max resid 0.07692795
Run 10 stress 0.1709846
Run 11 stress 0.1508712
Run 12 stress 0.1496642
... Procrustes: rmse 0.01285554 max resid 0.07716671
Run 13 stress 0.1498306
... Procrustes: rmse 0.05243516 max resid 0.3271674
Run 14 stress 0.1496489
... Procrustes: rmse 0.01266723 max resid 0.07694581
Run 15 stress 0.149831
... Procrustes: rmse 0.05247949 max resid 0.3275625
Run 16 stress 0.1573579
Run 17 stress 0.1496498
... Procrustes: rmse 0.01268426 max resid 0.07733981
Run 18 stress 0.1635248
Run 19 stress 0.1498309
... Procrustes: rmse 0.05245881 max resid 0.3273925
Run 20 stress 0.1496495
... Procrustes: rmse 0.01267903 max resid 0.07724868
*** Solution reached
jac_nmds$stress #Check the stress. Less than 0.1 is good. Less than 0.5 is better. This will be different each time, since it is iteratively finding a unique solution each time (although the should look similar)
[1] 0.1495161
# Extract points from nmds and merge into data frame with metadata
jac_nmds_df <- data.frame(jac_nmds$points) %>%
rownames_to_column(var = "SampleID")
# Merge metadata into the pcoa data table
jac_nmds_df <- left_join(jac_nmds_df, metadata, by = "SampleID")
head(jac_nmds_df)
## Plotting euclidean distance NMDS
nmds_jaccard <- ggplot(jac_nmds_df,aes(x = MDS1, y = MDS2, color = Station, shape = Bayside)) +
geom_point(size = 4) +
scale_color_brewer(palette="Paired") +
theme_bw() +
labs(x = "NMDS 1", y = "NMDS 2", title = paste0('Jaccard Distance NMDS, Stress = ', round(jac_nmds$stress,2))) +
coord_fixed(ratio = 1)
nmds_jaccard
ggsave("figures/nmds_jaccard.eps",nmds_jaccard, width = 7, height = 5, units = c("in"))
This is still a relatively high stress (>0.1) so should be interpreted with caution. Similar to Aitchison-distance nMDS but there is a little more separation of East vs West on NMDS 2 axis.
Next try a Bray-Curis NMDS, which will represent differences in presence/absence among samples and relative abundance, thus emphasizing impacts of highly abundant taxa.
bray_nmds <- metaMDS(bray_dmat, k=2, autotransform=FALSE) # Run the ordination. Distance matrix was already calculated above
Run 0 stress 0.1628608
Run 1 stress 0.1498312
... New best solution
... Procrustes: rmse 0.08567004 max resid 0.321179
Run 2 stress 0.1511991
Run 3 stress 0.1496848
... New best solution
... Procrustes: rmse 0.01159486 max resid 0.06870893
Run 4 stress 0.1496492
... New best solution
... Procrustes: rmse 0.05214599 max resid 0.3270977
Run 5 stress 0.1805538
Run 6 stress 0.1495165
... New best solution
... Procrustes: rmse 0.01281364 max resid 0.07805366
Run 7 stress 0.1511997
Run 8 stress 0.1496853
... Procrustes: rmse 0.05087238 max resid 0.3273574
Run 9 stress 0.166339
Run 10 stress 0.1496489
... Procrustes: rmse 0.01281149 max resid 0.07792418
Run 11 stress 0.1508707
Run 12 stress 0.1568975
Run 13 stress 0.1496846
... Procrustes: rmse 0.05087779 max resid 0.3278401
Run 14 stress 0.1625699
Run 15 stress 0.1496844
... Procrustes: rmse 0.05084287 max resid 0.327484
Run 16 stress 0.1937921
Run 17 stress 0.1710636
Run 18 stress 0.1508709
Run 19 stress 0.1648571
Run 20 stress 0.1627921
*** No convergence -- monoMDS stopping criteria:
20: stress ratio > sratmax
Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls
4: In readChar(file, size, TRUE) : truncating string with embedded nuls
5: In readChar(file, size, TRUE) : truncating string with embedded nuls
6: In readChar(file, size, TRUE) : truncating string with embedded nuls
7: In readChar(file, size, TRUE) : truncating string with embedded nuls
8: In readChar(file, size, TRUE) : truncating string with embedded nuls
9: In readChar(file, size, TRUE) : truncating string with embedded nuls
10: In readChar(file, size, TRUE) : truncating string with embedded nuls
bray_nmds$stress #Check the stress. Less than 0.1 is good. Less than 0.5 is better. This will be different each time, since it is iteratively finding a unique solution each time (although the should look similar)
[1] 0.1495165
# Extract points from nmds and merge into data frame with metadata
bray_nmds_df <- data.frame(bray_nmds$points) %>%
rownames_to_column(var = "SampleID")
# Merge metadata into the pcoa data table
bray_nmds_df <- left_join(bray_nmds_df, metadata, by = "SampleID")
head(bray_nmds_df)
## Plotting euclidean distance NMDS
nmds_bray <- ggplot(bray_nmds_df,aes(x = MDS1, y = MDS2, color = Station, shape = Bayside)) +
geom_point(size = 4) +
scale_color_brewer(palette="Paired") +
theme_bw() +
labs(x = "NMDS 1", y = "NMDS 2", title = paste0('Bray-Curtis Distance NMDS, Stress = ', round(bray_nmds$stress,2))) +
coord_fixed(ratio = 1)
nmds_bray
ggsave("figures/nmds_bray.eps",nmds_bray, width = 7, height = 5, units = c("in"))
Very similar to Jaccard results. High-ish stress (0.15)
–> CONTINUE HERE
XXX
I am interested in pulling out which taxa define the different in PA and FL.
The PhILR package has some interesting capabilities to do this using sparse logistic regressions.
Load glmnet and fit sparse logistic regression model (Note, if deriving a variable based on another with more than 2 factors, you need some extra lines. See tutorial linked above for details)
library(glmnet)
glmmod <- glmnet(ps.philr, sample_data(ps)$Size_fraction, alpha=1, family="binomial")
Identify the balances that define Total vs. Pfree. Use the default parameters from the tutorial for now.
top.coords <- as.matrix(coefficients(glmmod, s=0.2526))
top.coords <- rownames(top.coords)[which(top.coords != 0)]
(top.coords <- top.coords[2:length(top.coords)]) # remove the intercept as a coordinate
Visualize the balance values with respect to Total/ Pfree for each of these 3 balances
ps.philr.long <- convert_to_long(ps.philr, get_variable(ps, 'Size_fraction')) %>%
filter(coord %in% top.coords)
ggplot(ps.philr.long, aes(x=labels, y=value)) +
geom_boxplot(fill='lightgrey') +
facet_grid(.~coord, scales='free_x') +
xlab('Size fraction') + ylab('Balance Value') +
theme_bw()
There are 3 balances which represent the differences. Find out what those are:
tc.names <- sapply(top.coords, function(x) name.balance(tree2, tax2, x))
tc.names
These balance names represent the highest node on the tree, which can include many taxa below it (which is why, for example, n5 is comparing Bacteria and Bacteria- there are different taxa under each of these nodes).
Check full taxonomy of each of the 3 balances- First n5:
votesn5 <- name.balance(tree2, tax2, 'n5', return.votes = c('up', 'down'))
votesn5[[c('up.votes', 'phylum')]] # Numerator at Class Level
votesn5[[c('down.votes', 'phylum')]] # Denomenator
And visualize the position of the n5 node on a tree:
library(ggtree)
tc.nn <- name.to.nn(tree2, top.coords)
tc.colors <- c('#a6cee3')
p <- ggtree(tree2, layout='fan') +
geom_balance(node=tc.nn[1], fill=tc.colors[1], alpha=0.6)
p <- annotate_balance(tree2, 'n5', p=p, labels = c('n5+', 'n5-'),
offset.text=0.15, bar=FALSE)
p
The above shows the split between the numerator of the balance (n5+: Cyanobacteria, Firmicutes, Fusobacteriota) and the denomenator (n5-: Acidobacteriota, Actinobacteriota, Armatimonadota, Bacteroidota, Chloroflexi, Cloacimonadota, Crenarchaeota, Gemmatimonadota, Marinimicrobia (SAR406 clade), Nitrospinota, Nitrospirota, Patescibacteria, Planctomycetota, SAR324 clade(Marine group B) and Verrucomicrobiota).
I am not positive what is represented by the numerator and denomenator here (Total or Pfree), although I can guess that n5+ [numerator] is Total because a) that is the first factor that appears in ps$Size_fraction and b) there are less phyla (3 in n5+ vs. 15 in n5-) and it is a smaller chunk on the tree. I will check the relative abundance of some of these phyla, facet-wrapped by size fraction, to confirm which balances are in higher [relative] abundance in Total vs. Pfree:
First make a function that plots pyloseq data as bars but doesn’t include those black lines (from here )
#
my_plot_bar = function (physeq, x = "Sample", y = "Abundance", fill = NULL, title = NULL,
facet_grid = NULL) {
mdf = psmelt(physeq)
p = ggplot(mdf, aes_string(x = x, y = y, fill = fill))
p = p + geom_bar(stat = "identity", position = "stack")
p = p + theme(axis.text.x = element_text(angle = -90, hjust = 0))
if (!is.null(facet_grid)) {
p <- p + facet_grid(facet_grid)
}
if (!is.null(title)) {
p <- p + ggtitle(title)
}
return(p)
}
# Check cyanobacteria, which is in the numerator for balance n5
ps_cyano_ra <- subset_taxa(ps_ra, phylum == "Cyanobacteria")
ps_cyano_barplot <- my_plot_bar(ps_cyano_ra) +
facet_wrap(~Size_fraction, scales = "free_x")
ps_cyano_barplot
These seem more abundant in the Total size fraction, which would be consistent with numerator being the “Total” fraction.
# Check Firmicutes, which is in the numerator for balance n5
ps_firmi_ra <- subset_taxa(ps_ra, phylum == "Firmicutes")
ps_firmi_barplot <- my_plot_bar(ps_firmi_ra) +
facet_wrap(~Size_fraction, scales = "free_x")
ps_firmi_barplot
Firmicutes reach higher abundances in Total size fraction, but on average I wouldn’t necessarily say they are different between size fractions (but remember this plot is looking at relative abundance, not composition, like PhILR)
# Check Fusobacteriota, which is in the numerator for balance n5
ps_fuso_ra <- subset_taxa(ps_ra, phylum == "Fusobacteriota")
ps_fuso_barplot <- my_plot_bar(ps_fuso_ra) +
facet_wrap(~Size_fraction, scales = "free_x")
ps_fuso_barplot
These also seem to be in higher abundance in the total, so I am beginning to decipher that the numerator of the balance is “Total” and the denomenator is “Pfree”.
Check some from the denomenator (which according to my logic should be in higher abundance in Pfree)
# Check Bacteroidota, which is in the denomenator for balance n5
ps_bacteroid_ra <- subset_taxa(ps_ra, phylum == "Bacteroidota")
ps_bacteroid_barplot <- my_plot_bar(ps_bacteroid_ra) +
facet_wrap(~Size_fraction, scales = "free_x")
ps_bacteroid_barplot
These relative abundances seem similar in both size fractions, but perhaps a little higher on average in the Pfree
# Check Verrucomicrobiota, which is in the denomenator for balance n5
ps_verruco_ra <- subset_taxa(ps_ra, phylum == "Verrucomicrobiota")
ps_verruco_barplot <- my_plot_bar(ps_verruco_ra) +
facet_wrap(~Size_fraction, scales = "free_x")
ps_verruco_barplot
These also seem similar in both size fractions, not necessarily higher in Pfree. Could be because these are relative abdundances and PhILR is assessing composition (the ILR transformed abudnance data). Check one group that has only very low rel abundance overall…
# Check Nitrospinota, which is in the denomenator for balance n5 and has lower abundance overall, so might be more obvious
ps_nitrospino_ra <- subset_taxa(ps_ra, phylum == "Nitrospinota")
ps_nitrospino_barplot <- my_plot_bar(ps_nitrospino_ra) +
facet_wrap(~Size_fraction, scales = "free_x")
ps_nitrospino_barplot
This is definitely in higher relative abundance in the Total, not what I was expecting… This makes me a bit suspicious. Some of these mismatches in interpretation of PhILR results and the relative abundances can be explained by differences in transformed (compositional) data and relative abundance. But here, this balance is barely present in the Pfree fraction. But it is in the denomenator of this balance comparison, so I am not confident about this method yet.
Check full taxonomy of n29
votesn1114 <- name.balance(tree2, tax2, 'n29', return.votes = c('up', 'down'))
votesn1114[[c('up.votes', 'order')]] # Numerator at order Level
votesn1114[[c('down.votes', 'order')]] # Denomenator at order level
And visualize the position of the n29 node on a tree:
tc.nn <- name.to.nn(tree2, top.coords)
tc.colors <- c('#a6cee3')
p <- ggtree(tree2, layout='fan') +
geom_balance(node=tc.nn[2], fill=tc.colors[1], alpha=0.6)
p <- annotate_balance(tree2, 'n29', p=p, labels = c('n29+', 'n29-'),
offset.text=0.15, bar=FALSE)
p
This region is a very small part of the tree, within the Gammaproteobacteria. Hard to decipher from the full tree.
Check the relative abundance of some of the identified groups from balance n29 in Total vs. Pfree. First CCM19a, which is in the numerator:
ps_CCM19a_ra <- subset_taxa(ps_ra, order == "CCM19a")
ps_CCM19a_barplot <- my_plot_bar(ps_CCM19a_ra) +
facet_wrap(~Size_fraction, scales = "free_x")
ps_CCM19a_barplot
Definitely in higher abundance in Total and is present in numerator (n29+)
Check something from the denomenator, CCD24
ps_CCD24_ra <- subset_taxa(ps_ra, order == "CCD24")
ps_CCD24_barplot <- my_plot_bar(ps_CCD24_ra) +
facet_wrap(~Size_fraction, scales = "free_x")
ps_CCD24_barplot
So while the abundance of CCD24 is higher in Pfree than CCM19a is, CCD24 is still in higher relative abundance in Total than Pfree. I think there might be a relative abudnance bias against finding the Pfree- favored groups, that is being oversome by the PhILR method, which is why I keep seeing this. Try one more from the denomenator,EPR3968-O8a-Bc78
ps_EPR3968_ra <- subset_taxa(ps_ra, order == "EPR3968-O8a-Bc78")
ps_EPR3968_barplot <- my_plot_bar(ps_EPR3968_ra) +
facet_wrap(~Size_fraction, scales = "free_x")
ps_EPR3968_barplot
Here, the relative abundance definitely seems higher in Pfree, which is consistent with this group being in the denomenator.
ALDEx2 is a compositional approach to differential abundance that is similar to ANOVA. It models the probability of observing a count as a log-ratio transformed probability distribution rather than as direct counts, therefore correcting for differences in sequencing depth and being insensitive to the removal of ASVs from a sample.
library(ALDEx2)
The ALDEx t-test can be run in one step but is doing a lot of things. Here’s a quick summary from Nicholas Olberding’s website:
aldex2_ps <- ALDEx2::aldex(data.frame(phyloseq::otu_table(ps)), phyloseq::sample_data(ps)$Size_fraction, test="t", effect = TRUE, denom="iqlr")
Plot effect size
#Plot effect sizes
ALDEx2::aldex.plot(aldex2_ps, type="MW", test="wilcox", called.cex = 1, cutoff = 0.05)
Differentially abundant taxa are those where the difference exceeds dispersion. Points in the top are more abundant in the Total samples and towards the bottom are more abudnance in the Pfree samples. Highlighted in red are taxa with BH-FDR corrected p-values.
However, the authors of ANDEx2 recommend using effect size rather than p-values:
# Modify taxonomy table
ps_taxonomy <- tax_table(ps) %>%
as("matrix") %>%
tibble::as_tibble(rownames = "ASV")
#Clean up presentation of table
sig_aldex2 <- aldex2_ps %>%
rownames_to_column(var = "ASV") %>%
filter(wi.eBH < 0.05) %>%
arrange(effect, wi.eBH) %>%
dplyr::select(ASV, diff.btw, diff.win, effect, wi.ep, wi.eBH)
# Match ALDEx table to taxonmy table and view
sig_aldex2 <- left_join(sig_aldex2, ps_taxonomy, by = "ASV")
sig_aldex2
Table is arranged by diff.btw and effect. diff.btw is the median difference in clr between the Total and Pfree groups. Effect is the the diff.btw/ max(diff.win) for all instances. According to the authors, effect size is the key metric here: it is robust for a variety of distributions and tells us which ASVs are reproducibly different between groups. It is also
The taxa with negative effect values are those with higher abundance in Pfree and those with positive effect values are those in higher abundance in Total (based on plot above). There is some inconsistency here with the PhILR results-
Negative values (Pfree): Several ASVs from Bacteroidota, several from Gammaproteobacteria (Burkholderiales), 2 SAR11 clades, 1 Actinobacteria. Positive values (Total): there are many but here is a brief summary Several different Alpha and Gammaproteobacteria, including the EPR3968-O8a-Bc78 and Steroidobacterales ASVs that came up before, several Chloroplast ASVs, several Desulfobacterota, several Bacteroidota, several Chloroflexi, several Nitrospirota, several Verrucomicrobiota, and more…
While the chloroplasts seem to be consistenly in the Total size fraction for both the PhILR and ANDEx2 approaches (which makes sense, since these are eukaryotic), other groups like Nitrospirota, Verrucomicrobiota, Gammaproteobacteria etc. are in the Total group by ALDEx method but in the Pfree according to the PhILR method. I think the ALDEx2 method is more widely accepted for differential abundance analysis of microbiome data. Still, I should check at least one more approach to be confident…
Above step takes a bit (20-30 mins) so save variables again up to this point incase
Since that step takes a long time, save all variables up to this point in environmnet as RData object
# save.image("HREmicrobiome_analysis_vars_upto_aldex2.RData")
If variables are lost, come back and load everything up to now:
# load("/Volumes/MyPassport/HREMicrobiomeproject2019_LargeFiles/dada2_decipher_analysis_files/post_analysis_R/HREmicrobiome_analysis_vars_upto_aldex2.RData")
Note on the above: Copied the Rdata object into my external storage. Load it from there.
Following Mandal et al. 2015 and Kaul et al. 2017. ANCOM is an Aitchison log-ratio approach to differential abundance. ANCOM-II is the same approach as ANCOM, with a pre-processing step that corrects for excessive zeroes (which most microbiome data have). Here I am following the tutorial from here, specifically the moving_pics.R tutorial in the scripts directory.
Clone the ANCOM repository (in terminal)
git clone https://github.com/FrederickHuangLin/ANCOM.git
library(exactRankTests)
library(nlme)
source("ANCOM/scripts/ancom_v2.1.R")
Read in data and do some pre-processing for correct format. Based on Kaul et al. 2017, the types of zeroes we would have in this dataset should be considered structural zeroes.
feature_table <- as.data.frame(otu_table(ps)) # absolute abundance table; taxa as rows
meta_data <- as.data.frame(sample_data(ps)) # meta data table
colnames(meta_data)[1] <- "Sample.ID" # simplfy column name of sample IDs from meta_data
sample_var = "Sample.ID"; # indicate this name in a variable
group_var = NULL # required for detecting structural zeroes and outliers.
out_cut = 0.05; # observations with proportion of mixture distribution less than out_cut will be detected as outlier zeros
zero_cut = 0.90; # Taxa with proportion of zeroes greater than zero_cut are not included in the analysis
lib_cut = 1000; # Samples with library size less than lib_cut are not included in the analysis
neg_lb = TRUE; # TRUE indicates a taxon would be classified as a structural zero in the corresponding experimental group using its asymptotic lower bound. More specifically, neg_lb = TRUE indicates you are using both criteria stated in section 3.2 of ANCOM-II to detect structural zeros
Pre-processing- to detect structural zeroes
prepro = feature_table_pre_process(feature_table, meta_data, sample_var, group_var,
out_cut, zero_cut, lib_cut, neg_lb)
feature_table = prepro$feature_table # Preprocessed feature table
meta_data = prepro$meta_data # Preprocessed metadata
struc_zero = prepro$structure_zeros # Structural zero info
ANCOM step- takes a while so comment out for now
# main_var = "Size_fraction";
# p_adj_method = "BH"; # default method is Benjamini-Hochberg procedure
# alpha = 0.05 # significance cutoff
# adj_formula = NULL; rand_formula = NULL
# t_start = Sys.time()
# res = ANCOM(feature_table, meta_data, struc_zero, main_var, p_adj_method,
# alpha, adj_formula, rand_formula)
# t_end = Sys.time()
# t_run = t_end - t_start # around 30s
#
# write_csv(res$out, "ancom_results.csv")
t_run
Above step took >24 hours so save variables again up to this point incase
Since that step takes a long time, save all variables up to this point in environmnet as RData object
save.image("HREmicrobiome_analysis_vars_upto_ancom.RData")
If variables are lost, come back and load everything up to now:
# load("/Volumes/MyPassport/HREMicrobiomeproject2019_LargeFiles/dada2_decipher_analysis_files/post_analysis_R/HREmicrobiome_analysis_vars_upto_ancom.RData")
Note on the above: Copied the Rdata object into my external storage. Load it from there.
Plot- FIGURE THIS OUT
# Step 3: Volcano Plot
# Number of taxa except structural zeros
n_taxa = ifelse(is.null(struc_zero), nrow(feature_table), sum(apply(struc_zero, 1, sum) == 0))
# Cutoff values for declaring differentially abundant taxa
cut_off = c(0.9 * (n_taxa -1), 0.8 * (n_taxa -1), 0.7 * (n_taxa -1), 0.6 * (n_taxa -1))
names(cut_off) = c("detected_0.9", "detected_0.8", "detected_0.7", "detected_0.6")
# Annotation data
dat_ann = data.frame(x = min(res$fig$data$x), y = cut_off["detected_0.9"], label = "W[0.9]")
fig = res$fig +
geom_hline(yintercept = cut_off["detected_0.9"], linetype = "dashed") +
geom_text(data = dat_ann, aes(x = x, y = y, label = label),
size = 4, vjust = -0.5, hjust = 0, color = "orange", parse = TRUE)
fig
#ggsave("images/moving_pics.jpeg", height=5, width=6.25, units='in', dpi = 300)